Dealing with entanglement of continuous variables: 
Schmidt decomposition with discrete sets of orthogonal functions 
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We propose a method for obtaining the Schmidt decomposition of bipartite systems with contin- 
uous variables. It approximates the modes to the prescribed accuracy by well known orthogonal 
functions. We give some criteria for the control of errors. We illustrate the method comparing 
its results with the already published analysis for entanglement of biphotons. The agreement is 
excellent. 
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I. INTRODUCTION 

Bipartite and multipartite entanglement is one of the 
features that give rise to many of the developments of 
quantum computation and information, like quantum 
teleportation ll Q and quantum cryptography 0, 0, 
among others f| . The evaluation of the entanglement of 
a composite state is thus a main task to be done. For this 
purpose, the Schmidt decomposition 0,0 nas proven to 
be a valuable tool, for systems with just two components. 

In this paper we consider the case of continuous vari- 
ables entanglement. For us, these variables may be 
{a+at, i(aJ—a)} which commute as phase space variables 
do. We also refer to continuous variable entanglement in 
systems described by momentum and/or energy obscrv- 
ablcs. Precisely, the entanglement of continuous vari- 
ables stems from the original EPR article However, 
the treatment of systems with continuous variables is far 
from straightforward. Until now, Schmidt decomposition 
in the continuous case re quir ed solving the correspond- 
ing integral equations 0, E3, El E 03 They had to 
be discretized, losing the continuous dependence of the 
initial state. Here we propose a method to perform the 
Schmidt decomposition for this case, to the accuracy de- 
sired, keeping the continuous character of the variables. 
This method consists of two steps: 

1) We decompose the bipartite system wave function, 
f(p,q), by using two discrete and complete sets of or- 



thogonal functions, 
form: 



{oi 1} (p)}, {or(g)>, of L 



(2), 



f(p, q ) = Y,C m nO^(p)Oi?( q ) 



the 



(1) 



The purpose of this step is to transform the continu- 
ous problem into a discrete one (a necessary step for the 
numerical computation) , while preserving the continuous 
dependence of f(p, q). 
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2) Then we apply the (finite dimensional) Schmidt pro- 
cedure to 0} in order to write the wave function f(p, q) 
as diagonal sum of biorthogonal terms: 



(2) 



The orthogonal functions ipn\p), ipn\q) -the modes- 
will be some particular linear combinations of On\p), 

On {q) , respectively. Notice that we are using the 
Schmidt procedure for discrete systems to obtain the de- 
composition for the continuous case. This is much more 
tractable, as it implies diagonalizing matrices instead of 
solving integral equations. 

The rationale for this procedure is the expectation that 
only a few O n will suffice: A handful of appropriate or- 
thogonal functions will approximate f(p, q) to the desired 
accuracy. We finish by pointing out some properties of 
this method, namely 

a) We obtain complete analytic characterization of the 



modes ipn\p), iftn' (?) to the desired precision. Our 
method surpasses the standard numerical procedures in 
that keeps the continuous features present in f(p, q). 
b) We remark the portability of the attained modes 



:,(2) 



' l Pn\p)^ ^n'{q) that are ready for later uses. 

c) For the physical system analyzed in this paper 
(biphoton), we found that with 26 x 26 C mn matrices 
the error was of around 2%, that is, the number of O n 
functions required is small. For other systems studied 
that we do not include here, the convergence was even 
better: For the case of two electrons which interact elec- 
trostatically, the obtained error with 12 x 12 matrices was 
of 0.7% (Schmidt number K = 2.4). 

In this paper we begin with a detailed exposition of 
our method in Sect. [H] Then, in Sect. IIIII we apply it to 
a relevant case: the entanglement of two photons created 
by parametric down-conversion. We compare our results 
(leading to well known, continuous functions) with those 
computed by numerical methods [lfjj that produce sets of 
points: discrete functions. Both methods agree remark- 
ably well. Finally, in Sect. IIVI we decompose the Dirac 
delta, a case of physical and mathematical interest. 



.,(2) 
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II. SCHMIDT DECOMPOSITION WITH 
DISCRETE SETS OF ORTHOGONAL 
FUNCTIONS 

We consider a bipartite quantum system formed by 
two subsystems Si and S*2. Some examples are two pho- 
tons entangled by parametric down-conversion, a photon 
emitted by an excited atom and as a result entangled 
with it or two charged particles which interact electri- 
cally. This system is described by the vector state 



|0)= /dpdg/(p,g)o[ 1) (p)oj a) (?)|0,0) (3) 

(\\f(p,q)\\ 2 = [ dpdq\f{p,q)\ 2 < oo) 



where at^(p), a\ 2 ,(q) are the creation operators of a 
particle associated to the subsystems Si and S2 . p and q 
are continuous variables associated to Si and S2 respec- 
tively, which can represent momenta, energies, frequen- 
cies, or the like. In general, the analysis is made in an 
ad hoc kincmatical situation in which p and q turn out 
to be one-dimensional variables, p € (ai, &i), q £ (a 2 , b 2 ). 
In the following we assume this is the case. In addi- 
tion, there can be discrete variables (like the spin) to be 
treated with the Schmidt method, that we do not include 
here to avoid unwieldy notation. 

Our method works as follows: 

We consider two denumerable, complete sets of orthog- 
onal L 2 functions {On (p)}> {0„\q)} n = 0,1,. ..,00, 
each one associated to each particular subsystem S a 
(a = 1,2). These functions obey 



(4) 
(5) 



dkO£>(k)Oi a \k) = 5 rnn 
Y,0^*{k)0^{k') = 8{k-k') 



1) Our first step is to expand the wave function f(p, q) 

as a linear combination of the O n °^ , translating the con- 
tinuous problem into a discretized one. Thus we work 
with the discrete coefficients of the linear combination, 
though the continuous character of the state is preserved 



in the k dependence of the O 
reads: 



(a) 



functions. The expansion 



f(p,q)= £ C mn O%(p)0%\q) (6) 

m,n— 

where the coefficients C mn are given by 



C„„ = / 1 dpO$*(p) [' dqO£>(q)f(p,q) (7) 



2) Our second step is to apply the Schmidt decompo- 
sition to the discretized bipartite state ©, as is usually 
done for finite dimension Hilbert spaces (diagonalizing 
matrices, instead of solving integral equations). In order 
to do this, it is necessary to truncate the expansion JfjJ, 
something that is possible to a certain accuracy due to 
the fact that J dpdq\f(p,q)\ 2 < 00 (f(p,q) is in principle 
normalizable), and the expansion is in orthogonal func- 
tions, so the coefficients C mn go to with increasing to, n 
(see below). 

We truncate the series © at to — mo, n = no, with 
tno < no, without loss of generality. The Schmidt proce- 
dure leads to where 



•71 = 

mo no 



0, ...,m 



(8) 
(9) 



Here the matrix V is the (transposed) eigenvectors ma- 

tnxofMy =m; 4 = E"L C m C7„: 



mo 



Xj Vji 



(10) 



and {A;}i = o,...,m are the eigenvalues of M. 
There are two sources of error in this procedure: 

a) Truncation error: This is the largest source of er- 
ror in our method. Inescapably, the series 10 must end 
at some finite m, n when attempting to obtain some 
specific result. This step is possible to a certain ac- 
curacy because the function f(j>, q) is square- integrable 
and we are expanding it into orthogonal functions, so 
Em=oE,T=o |C mn | 2 < 00 and thus C mn -> when 

TO, 71 — * OO. 

The particular choice of the orthogonal functions 
will affect how fast the C mn go to zero. Hence, the elec- 
tion of these functions for a particular physical problem 
will be a delicate task. To reach the same accuracy with 
different sets {O^} it will be necessary in general to 
consider a different pair of cut-offs {too, no} for each of 
the sets. 

b) Numerical error: This is a better controlled source. 
It includes the error in calculating the coefficients C mn 
via and the one produced when diagonalizing the ma- 
trix M = CCl 

The suitable quantity to control the convergence for 
a particular f(p,q) and a specific set {O^} is the well 
known (squared) distance no between the function 
f(p, q) and the Schmidt decomposition obtained with 
cut-offs {mo, no} (mean square error): 



3 



^\p)^\q)\ 2 



(11) 



this expression gives the truncation error. It will go 
to zero with increasing cut-offs according to the specific 
{O^} chosen. 

Another easily computable, less precise way of control- 
ling the convergence is given by the fact that (with no 
cut-offs) Em=o A ™ = \\f(P^)\\ 2 and thus 



"m ,n 



1 - 



Em 
m— 



Wf(p,qW 



(12) 



is other measure of the truncation error, where here A m 
is calculated with cut-offs {mo, no}. Would we compute 
the A„ exactly, then d 1 — d 2 . In practice this can not 
be done because our A„ are the eigenvalues of the mo x 
too matrix My, that depend slightly on mo, no. Both 
distances behave in a very similar way, as we show in 
Fig. n an d Fig- HI though d 2 is more easily computable 
than d 1 . 

The election of the two sets of orthogonal functions 
for a particular physical problem, {O^ } Q =i,2 can be 
approached from two different points of view, according 
to the feature one desires to emphasize: fundamental or 
practical. 



A. Fundamental point of view 

The election of the orthogonal functions in a particular 
problem can be done according to the specific intervals in 
which the variables p, q take values for that case. Typical 
examples of discrete sets of orthogonal functions are the 
orthogonal polynomials, defined in a variety of intervals. 
For example, a possible choice to describe one dimen- 
sional momenta p e (—00,00), are the Hermite polyno- 
mials, On\p) ~ H n (p). The equivalence sign indicates 
here that the polynomial must be accompanied by the 
square root of the weight function in order to be cor- 
rectly orthonormalized, and normalization factors must 
be included. If, on the other hand, the variable of inter- 
est in a specific problem is bounded from below, like the 
energy of a free massless particle p € (0,oo) , then the 

election could be Laguerre polynomials, On\p) ~ L n (p). 

The criterion for choosing the orthogonal functions 
0( a ) according to the intervals in which p, q are defined 
has a fundamental character. For example, the localiz- 
ability in configuration space of the Fourier transforms of 
the modes (JHJ), ©, depends critically on the intervals in 
which these modes are defined [lj| . Only if we choose the 
functions to be defined exactly in the same intervals 
as the amplitude f(p,q), may the Fourier transforms of 
the modes have the right localization properties. In spite 



of that, this point of view may not be the most suitable 
one, as it may give slower convergence than the point of 
view presented below. 



B. Practical point of view 

In this case, the election is approached with the goal 
of improving the convergence. The are chosen 

here according to the functional form of f(p,q). The 
closer the lowest modes are to / the lesser the num- 
ber of them necessary to obtain the required accuracy. 
What we are looking for here are that maximize 

dpOg ] *{p) j£ dqO^*(q)f{p, q) for low to, n. 

In some cases this practical point of view will be more 
useful than the fundamental one. For example, suppose 
the amplitude for a particular problem is of the form 
f(p,q) = g(p,q)e~( p / rTp ' > e^^l *) j with g(p,q) a slowly 
varying function of p, q. In this particular case it is rea- 
sonable to choose the functions as Hermite poly- 
nomials, because their weight functions are indeed gaus- 
sians. This leads to On^ (p) oc H n (p)e~ p I 2 , and similarly 
for O i2) (q). 

We call this approach practical because the conver- 
gence is reached faster. There is a price: the information 
relevant to localization might be lost. 



III. ENTANGLEMENT OF CONTINUOUS 
VARIABLES IN PARAMETRIC 
DOWN-CONVERSION 

In this section we consider a realistic case of bipho- 
tons already studied in the literature [lfi : two pho- 
tons entangled in frequency through parametric down- 
conversion. We apply our method to this physical sys- 
tem in order to obtain the Schmidt decomposition and 
the structure of modes without losing the analytic char- 
acter within the target accuracy. 

The system under study is a biphoton state generated 
by parametric down-conversion (PDC) of an ultrashort 
pump pulse with type-II phase matching. The amplitude 
in this particular case takes the form |10| 



f(uj ,u e ) = exp[-(w + ui e - 2uj) 2 /a 2 } 
xsinc{i[(w - uj){k' -k) + (uj e - Ld)(k' e - k)]/2} (13) 

where uj ,oj e € (0, 00) are the frequencies associated 
to the ordinary and extraordinary fields respectively, k' a 
and k' e are the inverse of group velocities at the frequency 
Co, k is the inverse group velocity at the pump frequency 
2d), L is the PDC crystal length and a is the width of 
the initial pulse. Typical values for these parameters are 
{k-k' e )L = 0.213ps, (k-k' )L = 0.061ps, u = 2700ps-\ 
L = 0.8mm and a — 35ps _1 . 

We perform now the following change of variables 
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5 10 15 20 25 

m (= n ) 



FIG. 1: d i m(] na as a function of the cut-offs {mo, no}. 

P = L p = (k-k' )La (14) 

g = — ~ UJ ; L q = (k-k' e )La (15) 
a 

and thus obtain 

/(p, g) = e-^+tfsmcl(L p p + L q q)/2] (16) 

We have applied our method to the function l|16fl (once 
normalized) according to section^ in the following way: 

We choose as orthogonal functions Hermite polynomi- 
als, because their weights are gaussians and a gaussian 
appears in l|lb|) . These polynomials were used in (la ] for 
PDC in some particular cases which are exactly solvable. 
The orthonormal sets we chose, looking for maximizing 
the C mn for the lowest to, n, were 

0^\k) = {^2 n n\)-^ 2 H n (k)e-^l 2 a = 1,2 (17) 

This choice of polynomials is suitable for the practical 
approach (subsection ^ Bp . taking into account that u> 3> 
a and thus the interval of definition of f(uj 0l uj e ) can be 
restricted to a region centered in u> of width ~ a in u> , u> e . 
We did a careful analysis of this, that for brevity we do 
not show here. Notice that our conclusions would not 
apply in the case L p = L q . 

We have considered cut-offs too — no taking values {5— 
25} and followed the steps of Sect. [H] We have computed 
the eigenvalues X n of the Schmidt decomposition for 
each pair {mo, no}. We have also computed the modes 
© and ©• 

In Fig. ^ we plot the distance d mo „ (|11|) as a function 
of too = no, to show how fast the convergence is. With 
mo = no = 25 the truncation error is of 2%. We also 
plot in Fig. |21 the distance ^ o no which serves 

as another measure of the convergence, as a function of 
rn = n . We obtained d\ 5 25 = 2%. 

Regarding now the most precise case considered, mo = 
no = 25, we plot in Fig. [3]the eigenvalues \ n for different 
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FIG. 2: djL, „„ as a function of the cut-offs {mo, no}. 
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FIG. 3: Eigenvalues \ n versus index n. 



values of n, observing good agreement with the results 
existing in the literature [l(j ■ For this case we also plot in 
Fig. @]the modes JSJ) and 10 for i = 0, 1, 2, 3, confirming 
the validity of the method when comparing with [lfj . 
The modes are given explicitly by: 

25 

^(fc) = e- k2 /*J2(^2 n nr 1/2 A%lH n (k) (18) 

71=0 

m = 0, ...,25 a = 1,2 

where the values of the coefficients Ami are obtained 
through JSJ and JjJJ • The actual properties of the modes 
(|18fl depend on these values. In fact, the parity and num- 
ber of nodes is determined by them, taking into account 
that H n is a polynomial of degree n, parity (— ) n and 
having n nodes. 

A good approximation to the ipo (p) obtained with our 
procedure is 

tp^ (p) = e- p2/2 (0.81395 - 0.14764p 2 + 0.00821p 4 ) (19) 
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FIG. 4: Modes ip£ \p), i/>i 2) (g) as a function of p = and 
g= for n = 0,1,2,3. 



m = n 


(3 = 1.0 


13 = 0.5 


(3 = 2.0 


25 


0.020 


0.13 


0.037 


20 


0.024 


0.19 


0.041 


15 


0.032 


0.27 


0.050 


10 


0.062 


0.38 


0.064 



TABLE I: t4 ,n for i 3 = 1-0,0.5,2.0 and m = n = 
25,20,15,10. 



0^(fc) = ^=^=H„(/3fc)e- (W2/2 a = 1,2 (21) 

We applied our method to the amplitude l|16|) with 
these sets of orthogonal functions, for (3 = 1.0,0.5,2.0, 
and cut-offs m = n = 25, 20, 15, 10. We show in Table 
H|the values of d 2 ngno for these specific parameters. 

Clearly, the convergence is better for the case (5 = 1, 
which we used in the preceding calculations. In case we 
chose another type of orthogonal function for (|16(l (La- 
guerre, Legendre,...), the convergence would have been 
much worse because of the specific shape of that ampli- 
tude. 



This expression has a deviation (squared distance) of 
10~ 5 from the whole mode obtained including terms until 
p 25 , which is the greatest power appearing for mo = uq = 
25. On the other hand, d\ A - 4 5)2 5 = 0.213 > 10~ 5 . 
From Q19[l it can be seen that in this mode the even com- 
ponents are greater than the odd ones (these are negligi- 
ble), so it is an even state, as shown in Fig. 0] 

Another example is the approximation to if)\ (q) 

%pf\q) = e - 92/2 (2.91088g - 3.54070g 3 + 1.29062g 5 
-0.20402g 7 + 0.01598g 9 - 0.00063g n + O.OOOOlg 13 ) 

(20) 

This has a deviation (squared distance) of 10~ 4 from 
the whole mode obtained including terms until q 25 . On 
the other hand, d} 3jl3 - d\ 5 25 = 0.020 > 10~ 4 . More 
terms are needed in (|27)|) . because they go to zero more 
slowly with increasing powers of q. Here the most im- 
portant components are the odd ones (the even ones are 
negligible), leading to an odd parity state, as shown in 
Fig.EI 

To show how the convergence of the method de- 
pends on the specific family pairs of orthogonal functions 

{On^ (p)}, {On^ (q)} chosen, we consider the cases of Her- 
mite orthogonal functions depending on a parameter j3 
related to the width of the gaussian, fixed for each family 
pair: 



IV. MAXIMUM ENTANGLEMENT: THE 
DIRAC DELTA 

Another interesting case is the Dirac delta. Here we 
have f(j),q) — 8{p — q) and we take the same interval 
(a, b) for p and q. We consider complete sets of orthonor- 
mal functions satisfying On (k) — On^*(k). A particular 
case is when they are real functions, as for example the 
typical orthogonal polynomials (Legendre, Hermite, La- 
guerre, Chebyshev,...) are. We must take into account 
that the Dirac delta is not a function but a distribution, 
and indeed is outside L 2 . However, we can calculate the 
C mn and study how much entanglement does this state 
have. We obtain straightforwardly C mn = S mn . This 
gives 



5(p-q) = Y, n h ^ nH<l) (22) 

n=0 

which is just the resolution of the identity as given in 
(J5J. The Schmidt decomposition of the Dirac delta is 
not unique, because all the weights vA™ are equal to one 
(they are degenerate). In fact, the decomposition can be 
done with any complete, denumerable set of orthonormal 
functions, in the form (J5J. This expression can be seen 
as an infinite entanglement case, in the sense explained 
below. The fact that all the weights are equal to one, 
makes sense only because we are considering a distribu- 
tion, not an L 2 state. The sum of the squares of the 
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weights, which must be equal to the square of the norm 
of the function f(p, q), diverges because the Dirac delta 
is not square-integrable. 

A possible measure of the entanglement of a state 
f(jp, q) in its Schmidt decomposition J3J) is given by the 
von Neumann entropy |a] 



S = - ^2 A„ log 2 A„ 



(23) 



This is usually called the entropy of entanglement. 

The state of L 2 closer to (122(1 is the case of an entangled 
state with N diagonal terms with equal A„, when N goes 
to infinity. To be correctly normalized it verifies A„ = 
1/N, n = 0, ...,N-1 and 



Urn 

JV->o 



N-l 1 1 

n=0 



lim log 9 — 



00 



(24) 

This is the maximum entanglement case. This provides 
an estimate of the entropy of the Dirac delta (were it in 
L 2 ). 

The origin of our interest in the 5 comes from the fact 
that f(p, q) may evolve in time towards a Dirac delta, 
as it happens in time dependent perturbation theory of 
quantum mechanics. The entanglement in these cases 
would increase with time towards its maximum, corre- 
sponding to the Dirac delta. 



V. CONCLUSIONS 



In this paper we have introduced a method for com- 
puting the Schmidt decomposition of a bipartite state 
with continuous degree s of freedom. In the existing liter- 
ature 0, 0, 0, 1 121 ll.il] the decomposition produced sets 
of points as approximation to the modes. Our method 
gives linear combinations of the well known orthogonal 
functions as approximation to them. When these func- 
tions are chosen properly, a handful of them is enough to 
reach the desired accuracy. We introduce some criteria 
for the control of convergence and truncation error. The 
result of our method for the decomposition of a bipho- 
ton state produced by parametric down-conversion agrees 
with the numeric results in the literature [lOj. We also 
touch on the last stage of evolution of entanglement for 
determined systems, analyzing the Dirac delta case. 
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